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A novel algorithm for estimating quark propagation for lattice hadron spectroscopy calculations 
will be presented. Lattice QCD is a framework in which the nonperturbative, first principles com- 
putation of hadronic correlation functions is possible by numerical simulations. It requires a lattice 
regulation of QCD in a finite Euclidean space-time where correlation functions are evaluated nu- 
merically via Monte-Carlo importance sampling methods. The asymptotic behaviour of correlation 
functions in Euclidan time is then used to extract the energies of the hadronic states of interest. 

This approach of computing correlation functions becomes numerically very challenging as one 
attempts to evaluate multi-hadron correlators and disconnected diagrams due to the rapid increase 
in the number of quark propagators involved in the calculation. Each quark propagator is obtained 
from an inversion of the 4-dimensional Dirac matrix in a finite, but large volume. Disconnected 
diagrams (quark loop diagrams) are particularly problematic because they require quark propagators 
from every point on a timeslice to every other point on the same timeslice, for all timeslices of the 
lattice. 

The cost of inverting the Dirac matrix, in terms of CPU cycles, increases as the quark mass being 
simulated approaches the physical values and the space-time volume is enlarged to reduce the finite 
size effects. It is not practical to compute these "all-to-all" quark propagators with the current 
resources available and this limits the physics that can be addressed in Lattice QCD. 

The standard solution is to stochastically estimate the all-to-all propagators with random noise 
sources. This method, however, introduces a lot of noise into the calculation which is reduced by 
performing more and more inversions. The stochastic LapH algorithm for inverting quark propa- 
gators avoids "using noise to cancel noise" approach by first reducing the space by cutting out the 
high frequencies (LapH) and then introducing noise sources in the LapH subspace. 

The usual volume scaling problems of all-to-all algorithms are absent in this method and the 
judicious choice of noise-dilution schemes makes the algorithm practical for real simulations that 
are performed near physical values of the quark masses. The signal for selected correlation functions 
will be shown to demonstrate the efficacy of the algorithm. 



I. INTRODUCTION 

There has been significant progress in lattice QCD simulations in recent years now that dynamical simulations 
can even be done at physical quark masses ([3], [Tj). Calculations of level splittings due to electromagnetic 
corrections (see review [5] and references therein) have also been attempted in recent years, as well as realistic 
attempts to calculate binding energies of atoms ([!]). The computation of the excited state spectra of baryons 
and mesons has been ongoing for the single particle sector in both Nf = 2 and Nf = 2 + 1 dynamical calculations 
at relatively large quark masses. The Hadron Spectrum Collaboration has been developing the tools necessary 
to complete this analysis starting with the design of the lattice interpolating operators for single particle states 

(El)- 

The question of mixings of excited single-particle states with multi-particle states has been pointed out in 
selected channels in recent spectroscopy papers (see for example [7], [T], [Sj and references therein), but these 
studies do not see multi-particle states in their spectrum even where they are expected (given the parameters 
of the simulations) . The main problem appears to be similar to the situation of the absence of string-breaking 
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of the static potential in dynamical simulations without explicit two-meson operators in the variational basis. 

In these proceedings, we report on the progress made in filling in this void by using an estimate of all-to-all 
quark propagators via the stochastic Lap _ff algorithm introduced earlier this year (|3]). We apply the algorithm 
to simple two-particle spectra with finite momenta hadron operators and to correlation functions which have 
disconnected diagram contributions, such as those in the isosinglet channel, including diagrams of mixing with 
scalar glueball operators. 

II. SPECTROSCOPY IN LATTICE QCD 

Spectroscopy in lattice QCD is conventionally carried out by computing hadron correlation functions and 
fitting the exponential decay of the signal in Euclidean time: 

(0|O(t)O(0)|0) ^ Aoe-^"* + Aie-^i* + • • • . (1) 

The interpolating operator O only needs to have a finite overlap with the state in question in the infinite time 
{t — ^ cxi) limit for the extraction of the ground state. An example of such an operator would be the single pion 
operator for which O — ip"/5ip is usually sufficient for extracting the pion mass. However, signal-to-noise ratio 
deteriorates rapidly for heavier states with increasing time so that interpolating operators with better overlaps 
are preferred to determine the energies with reasonable errors. 

The hadronic correlation functions are computed by contracting the appropriate combination of quark prop- 
agators to project out the required lattice irreps of the states being sought. It is the computation of these 
quark propagators (which requires the inversions of a large, four-dimensional matrix) that is the most time- 
consuming part of the analysis at light quark masses. It is usually not feasible to compute all elements of the 
quark propagator at the quark masses that we are interested in, even with the largest computational resources 
available to us today. Therefore, it is usually the case that only a handful of quark propagators are computed 
for spectroscopy calculations. 

The excited states in a given symmetry channel can be obtained by enlarging the interpolating operator basis 
and using the variational method to determine the optimal combination of those operators that project out the 
higher lying states. A correlation matrix is required in this case and effective interpolating operators must be 
included in the basis for this procedure to produce results with reasonable errors. For a two-particle state, it 
is preferred or even necessary that explicit two-hadron operators with finite momenta are included the basis of 
operators. 

Single-particle states with disconnected contributions in the correlation functions have the same problem as 
multi-particle states as quark propagators from many points on the lattice have to be computed to reduce the 
statistical errors that are typically intrinsically larger than connected correlation functions. The main problem 
comes from so-called "t-to-t" diagrams where a quark propagates out from timeslice t and propagates back to t 
again. The efficient evaluation of these diagrams is a necessary step towards the determination of the low-lying 
hadron spectrum from lattice QCD. The stochastic LapH algorithm was mainly designed to solve this problem. 
Details of the stochastic LapH algorithm are given in [3]. We outline the basic idea and its implementation for 
hadron spectroscopy calculations in the following sections. 

III. THE STOCHASTIC LAPH ALGORITHM 

The stochastic LapH algorithm is based on two main ideas that have been found to be effective in reduc- 
ing noise in the single timeslice-to-all and all-to-all propagators. The first is the idea of propagating only the 
low-lying modes of the quarks, as they are the modes which dominate the hadron correlation functions. These 
quark propagators go by the name of Lap _ff propagators (or distilled propagators, [9 ). LapH propagators are a 
particular implementation of the distillation algorithm with a ifeaviside step function cutoff of the eigenmodes 
{v^^\x,t)) of the the Laplacian operator. The main object which takes the place of the smeared quark prop- 
agators (called perambulators in [3]) take the following form in this case: V'''(t)]V[~^(t, to)V(to) where V(t) 

is the matrix formed from the column of eigenvectors of the Laplacian on timeslice t and M~^^ propagates the 
eigenmodes on timeslice to to t (and spin component (3 to a). We note that the perambulators enables one 
to inject finite momenta into the quarks but still connects timeslice to to t in the same way that the usual 
point-to-all propagators do, which can complicate the formation of hadronic correlation matrices. Moreover, 
the evaluation of t-to-t diagrams would require the computation of the perambulators from many timeslices on 
the lattice. 
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The second ingredient in the algorithm is the stochastic estimation of the quark propagators from all of the 
timeslices to make this a true all-to-all algorithm without computing all of the perambulators on the lattice. 
The use of stochastic noise to reduce the computational cost of all-to-all propagators is nothing new as it has 
been used in the computation of condensates and other quantities that have disconnected diagrams. In the 
stochastic LapH algorithm, .^4 noise is introduced in the LapH subspace in order to keep the noise out of the 
eigenvectors that were computed "exactly" . This means that each eigenvector of the Laplacian that we are 
using acquires a single noise factor instead of the usual A^space x ^coior x ^spin factors, greatly reducing the 
unwanted effects of the stochastic method. The noise vectors introduced in this way are denoted by (x, t) , 
where a denotes the eigenvector index, k denotes the spin index and r denotes the time index. The solutions 
of the Dirac equation with these noise sources carry similar indices, Lp^^^'^^x^t). 

The final ingredient in this algorithm is the use of dilution on the noise vectors for each quark line instead 
of averaging over a large number of independent sources to cancel the noise (see review in [lOj and references 
therein). This approach is much more effective than using multiple noise sources if one can find the right dilution 
scheme since the hopping of noise from one source to another is effectively removed with a judicious choice of the 
dilution scheme. Procedurally, the projection operators P^"*' are applied to the Z4 noise sources which dilutes 
the noise according to some scheme before they are used as the quark source in the Dirac equation. The diluted 
noise sources are denoted by, Q^'^\x,t) = P^'^"> g{x,t), where d is the dilution index and a similar index is used 
for the corresponding solution. 

The estimate of the full quark propagator for a given quark-line can then be written in the following compact 
form: 

d 

where the indices u and v are collective indices to simplify the notation and r denotes the particular Z4 noise 
source. A pion propagator, for example, can be formed from these propagators by performing the following 
contractions: 

a(t,io) =Tr[gV5Q'^75] . 
Similarly, a disconnected contribution to a meson correlation function can be written as 

C(i,i)=Tr[gL] . 

IV. PRELIMINARY RESULTS 

The preliminary results shown in this section were obtained on 2-\-\ dynamical, anisotropic lattices (anisotropy 
O'sjO't = 3.5) generated by the Hadron Spectrum Collaboration, Two values of the quark mass were used 
corresponding to pion masses of ~ 390 MeV (390 lattice) and ^ 240 MeV (240 lattice) with a spatial lattice 
spacing of ~ 0.12 fm. Several lattice sizes were generated with most of our preliminary results from the 16"^ x 128, 
20'^ X 128, 24'^ X 128 lattices and the number of configurations ranging from 100 ~ 300. 

A. J = 2 

The 1 = 2 sector is an exotic sector without any disconnected diagrams or box diagrams. Its simulation, 
therefore, has the longest history as all of the diagrams can be evaluated using the traditional point-to-all quark 
propagators. The main advantage of using the stochastic LapH method in this case is to increase statistics on 
a given lattice and to include finite momentum operators in the variational basis. An example of the effective 
mass of this state is shown in Figure [T] 

B. 1 = 1 

In the 7=1 sector, one must consider the mixing between the (single particle) rho and the two-pion state 
with sufficiently light quarks and large spatial volumes. The two pion diagram contains the t-to-t term as one 
of the quarks are exchanged between the two final state pions. The correlation matrix with a two-pion operator 
and a single rho operator must be formed to determine the mixing strength and the mass of the ground state, 
the An example for the signal of this state is shown in Figure [2] 
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FIG. 1: A sliding window plot of aoAI-^ for the 240 data plotted together with the negative of the effective mass (to show 
the quality of the data in the fit region). The dashed boxes are fits where the Xdof s-id Q did not satisfy the condition, 
Xdof < 1 ^-nd Q > 0-6- The two-loop xpt value is shown with a solid horizontal line and the dashed horizontal line is 
an estimate of the error due to the lack of the measurement of /jr on these lattices. 
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FIG. 2: The off-diagonal correlator of the p — >■ nn mixing term. 



C. 1 = 



The / = pion-pion correlation function contain both disconnected diagrams and box diagrams (which have 
t-to-t propagation on two of the "end timeslices"). The disconnected diagrams have much smaller magnitude 
than the connected diagrams (in this case) and are usually noisier than the connected diagrams computed with 
the same number of configurations. The problem in this channel is even worse as the complete study of the 



Proceedings of the DPF-2011 Conference, Providence, RI, August 8-13, 2011 



5 



isosingiet channel requires that we use a scalar glueball operator, quark-antiquark isoscalar meson-like operators, 
two-pion scattering operators and then study the mixing of the various components in this channel. The scalar 
glueball is one of the noisiest states (in the pure Yang-Mills theory) and the reliable extraction of the mixing 
with quark states has been a difficult problem to solve for some time. Since the stochastic LapH algorithm is a 
quark smearing algorithm, one may not immediately expect to gain anything in the glueball channel. However, 
the trace of the LapH operator has all of the symmetries needed for the 0"'"'" channel and can in principal be 
used as a glueball operator. Figure [3] shows the comparison of the effective masses between the usual single 
plaquette glueball operator and the LapH glueball operator. One sees that the LapH operator is just as good 
as the plaquette operator (or even better) for the scalar glueball. This implies that we can study all three types 
of / = operators with the stochastic LapH algorithm. The signal-to-noise is such that one can diagonalize 
this matrix to obtain the low- lying spectrum in the 0"*""'' channel (Fig. |4]). 
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FIG. 3: Comparison of the standard plaquette operator and the LapH operator for the scalar glueball. 
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FIG. 4: The optimized, diagonal correlators in the 1 = sector using a glueball, meson and multi-particle operators as 
the basis of operators in the variational basis. 



V. SUMMARY 



A new quark smearing algorithm for lattice spectroscopy was presented with preliminary results in the isospin 
0, 1 and 2 channels. A simple single exponential fit of the correlator of the state projected with two p — Q pion 
operators at M^r — 240 MeV gave results which were consistent with continuum chiral perturbation theory, [12] . 
Scattering states with finite momenta have also been obtained in the 1 = 2 channel at pion masses close to 
390 MeV and work is under way to compute the phase shift with lighter pions and larger volumes. This will 
enable us to directly compare phase shifts as well as threshold parameters. We have calculated mixing diagrams 
between the single-particle p operator and the two-pion scattering operators which show levels of noise that 
gives us confidence that decays can be studied with the stochastic LapH algorithm despite the presence of a 
t-to-t propagator. The most surprising result, however, may be in the isoscalar sector where the trace of the 
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LapH operator have allowed us to compute the mixing of a giueball operator, a scalar meson operator and two- 
pion operators with relatively low statistics. The preliminary results presented here on relatively small lattices 
strongly suggests that similar quality can be obtained at larger lattices and even lighter pion masses. We are 
currently generating the propagators needed to determine the low-lying spectrum of baryons and mesons with 
the correct treatment of the states above thresholds. 
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